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Abstract — The paper builds upon a recent approach to find 
the approximate bounds of a real function using Polynomial 
Chaos expansions. Given a function of random variables with 
compact support probability distributions, the intuition is to 
quantify the uncertainty in the response using Polynomial 
Chaos expansion and discard all the information provided about 
the randomness of the output and extract only the bounds 
of its compact support. To solve for the bounding range of 
polynomials, we transform the Polynomial Chaos expansion 
in the Bernstein form, and use the range enclosure property 
of Bernstein polynomials to find the minimum and maximum 
value of the response. This procedure is used to propagate 
Dempster-Shafer structures on closed intervals through non- 
linear functions and it is applied on an algebraic challenge 
problem. 

I. INTRODUCTION 

In interval analysis a fundamental problem is finding the 
interval bounds for the range of a real function. When such 
a function is monotone or it can be expressed in terms of 
arithmetic operations, then interval computations can be used 
to approximate the bounds of the response. However these 
bounds are gross overestimations due to the dependency and 
the wrapping effect [10]. 

Applications of interval methods can be found in estima- 
tion, optimization techniques, robust control, robotics and 
finance [15], [10], just to name a few. First, introduced by 
Moore [14], as a method to control for numerical errors 
in computers, the field of interval analysis has evolved 
with better approximations to the range of real functions as 
presented in Ref.[19], [18]. 

In the present paper we are interested in using interval 
methods to propagate epistemic uncertainty through nonlin- 
ear functions. In contrast to the aleatory uncertainty [26], 
[27], [11], defined by variability which is irreducible, the 
epistemic uncertainty is derived from incomplete knowledge 
or ignorance and can be reduced with an increase in in- 
formation. Due to their major differences, it is of great 
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importance that the two type of uncertainties to be modeled 
and propagated separately [6]. 

Unlike the probability theory, where the probability mass 
is assigned to singletons, in the Shafer's theory of evidence 
[23], the probability mass is assigned to sets, given its power 
on modeling ignorance. In conjunction with the Dempster's 
rule of combinations [4] which is a generalization of the 
of the Bayes' rule, the Dempster-Shafer (DS) theory of 
evidence offers a powerful methodology for representing and 
aggregating epistemic uncertainties. 

One can define Dempster-Shafer structures on focal ele- 
ments that are closed intervals on the real line for example. 
Person [7] shows how this structures can be transformed 
into probability bounds and vice-versa by discretization. To 
propagate these focal elements through system functions, 
it involves finding the solution to the interval propagation 
problem. 

We build upon a recent approach to find the approximate 
bounds of a real function using Polynomial Chaos expansions 
introduced by Monti [13], [24] and applied for worst-case 
analysis and robust stability. Given a function of random 
variables with compact support probability distributions, the 
intuition is to quantify the uncertainty in the response using 
Polynomial Chaos expansion and discard all the information 
provided about the randomness of the output and extract only 
the bounds of its compact support. 

Introduced by Norbert Weiner [28], Polynomial Chaos 
initially coined as the Homogeneous Chaos was used to rep- 
resent a Gaussian process as a series of Hermite polynomials. 
This method has been generalized to the Askey-scheme 
of orthogonal polynomials used to model random variables 
characterized by different probability density functions, in- 
cluding Beta and Uniform which have compact support [29]. 

The Polynomial Chaos is mathematically attractive due to 
the functional representations of the stochastic variables. It 
separates the deterministic part in the polynomial coefficients 
and the stochastic part in the orthogonal polynomial basis. 
This becomes particularly useful in characterizing the uncer- 
tainty of the response of a dynamical system represented by 
ordinary differential equations with uncertain parameters. 

To solve for the bounding range of polynomials, we 
propose to transform the Polynomial Chaos expansion in 
the Bernstein form, and use the range enclosure property of 
Bernstein polynomials to find the minimum and maximum 
value of the response [2]. The transformation does not 
require polynomial evaluations and it guarantees the global 
optimality of the bounds [8] and it is shown to be more 
efficient than existent interval global optimizers [20]. 



To demonstrate this approach for propagating Dempster- 
Shafer structures on closed intervals through nonlinear func- 
tions, we apply the proposed method on an algebraic chal- 
lenge problem used to investigate the propagation of epis- 
temic uncertainty [17]. 

The problem of propagating DS structures on closed 
intervals through nonlinear functions is stated in Section |II] 
and the proposed method is presented in Section |III] The 
numerical example is given in Section|IV]and the conclusions 
and future work are discussed in Section |V] 

II. PROBLEM STATEMENT 

A. Theory of Evidence 

The primitive function in the theory of evidence is the 
basic probability assignment (bpa), represented here by m, 
which is similar to the probability in the probabiUty theory. 
The bpa for a given set can be understood as the weight of 
evidence that the truth is in that set, evidence, which cannot 
be further subdivided among the members of the set. The bpa 
defines a map of the power set over \ht frame of discernment 
O to the interval [0, 1]: to : 2^ -> [0, 1]. The /oca/ element 
of m is every subset A C 51 such that m{A) > and the 
belief structure m verifies: 



Acn 



m{A) =■ 1 where m{Ai) = pi 



(1) 



In this work we are considering normalized belief struc- 
tures (closed-world assumption) which satisfy the following 
relation: m((/)) = 0, where (p is the null set. As an example 
consider the following body of evidence {Vt,m): Q, = 
{M,N,P} with m{{M,N}) = 0.3 and m{{N,P}) = 0.7. 

Based on the mass function two new functions can be 
induced. The belief function or the lower bound, Bel, which 
quantifies the total amount of support given to the set of 
interest A: 



Bel{A) = Y^ m{B) 



(2) 



BCA 



The plausibility function or the upper bound, PI, which 
quantifies the maximum amount of potential given to the set 
of interest A (here A is the complement of A): 

Pl{A) = Y^ m{B) = 1 - Bel{A) (3) 

The precise probability is bounded by the two quantities 
defined above, and when the equality is satisfied then the 
belief measure is just a probability measure and all the focal 
elements are singletons. 



Bel{A) <prob{A) < Pl{A) 



(4) 



Given two bpa's mi and 7TI2 based on independent argu- 
ments on the same frame of discernment, the Dempster's 
rule of combination provides the means to calculate the 
aggregation of the two belief structures: 

mi2iA^(t>) = T-^- y. mi{B)m2iC) (5) 



mi2{(t)) 



1- K ^ 

BnC=A 





where K = X]Bnc=<i "^i(^)"^2(C') represents the amount 
of probability mass due to conflict. 

While a number of combination rules have been derived 
to aggregate information [22], we present also the mixing 
rule which is used in the numerical example. Given n belief 
structures to be aggregated, the formula for the mixing rule 
is given by: 



I " 

TOl...„(A) = -Y WilTlilA) 

II -*- — ^ 



(6) 



1=1 



where wi are the corresponding weights proportional with 
the reliability of the sources. 

B. DS structures on closed intervals 

Given two nondecreasing functions F and F_, where 
F,F : R -)■ [0,1] and F{x) < F{x) for afl a; G R, we 
can represent the imprecision in the cumulative distribution 
function (CDF), F{x) ^ Prob{X < x), by the probability 
box (p-box) [F,F] as follows: F(x) < F{x) < F(x) [7]. 

A Dempster-Shafer structure on closed intervals can in- 
duce a unique p-box, while the inverse is not uniquely 
determined. Many Dempster-Shafer structures exist for 
the same p-box. Given the following body of evidence, 

{ {[Xi,Xi],pi) , {[X2,X2],P2) , •■• {[x^,Xn],Pn) }, the 

cumulative belief function (CBF) and the cumulative plausi- 
bility function (CPF) are defined by: 



CBF{x) = F{x) =Yp^ 

Xi<.X 

CPF{x) = F{x) =J2p^ 

x-<x 



(7) 
(8) 



Similarly one can obtain the complementary cumulative 
belief function (CCBF) and the complementary cumulative 
plausibility function (CCPF). 



CCBF{x) = 1 - CPF{x) = Y 



Pi 



CCPF{x) = 1 - CBF{x) = Y Pi 

Xi>X 



(9) 
(10) 



Thus the complementary cumulative distribution function 
(CCDF), Fc{x) = Prob{X > x), is bounded as follows: 



CCBF{x) < Fc{x) < CCPF{x) 



(11) 



Example: Consider the following body of evidence 
{ ([1,4], 2/3) , ([3,6], 1/3) }, the lower and the upper 
cumulative functions are plotted in Fig[T] 

C. Mapping of DS structures 

Consider the following function: y = /(a, b), where 
/ : R^ — > K, and a and 6 are given by the following 
bodies of evidence { {\a2_,ail_Pi) , •■• ([^,a^],K) } 
and { ([61, 6i],pJ) , ... {[bn, fo„],p^) } respectively. We are 
interested in finding the induced Dempster-Shafer structure 
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Fig. 1. P-box induced by the Dempster-Shafer structure 



in the y variable. The basic probability assignment describing 
y is given by [30]: 



mf{Y) 



E 



mi{Ai)m2{Bj) 



(12) 



where Y ^[y,y],Ai^ [o^, a^] and Bj 



[bj,bj 



Thus the problem of finding the mapping of a body of 
evidence on closed intervals is reduced to interval prop- 
agation [12]. This problem can be solved using the ad- 
vanced techniques developed in the interval analysis field 
[10]. However, due to the dependency problem the obtained 
bounds are conservative which is detrimental to the belief 
structure, since the evidence is assigned automatically to 
other elements which are not in the body of evidence. This 
problem becomes more acute when the uncertainty has to be 
propagated over a period of time. 

III. PROPOSED APPROACH 

We propose a new approach in approximating the prop- 
agation of intervals using a non-intrusive polynomial chaos 
method [3] in combination with Bernstein polynomials [25]. 
The approach of using polynomial chaos in propagating 
epistemic uncertainty has been considered previously in 
Ref.[13] for a limited number of arithmetic operations and 
relies on sampling or global optimization in the general 
case which is inaccurate for small number of samples and 
computationally expensive. 

A. Non-intrusive polynomial chaos 

The problem in solving the mapping of DS structures as 
given in Eq.(fT2b is to find Y = [y,y] such that: 



Y = fiA,B) 



(13) 



where A ~ [a, a] and B = [b,b]. While in this paper we 
present only the bivariate case, the method can be scaled up 
to the desired number of variables. 

The problem can be transformed into finding the stochastic 
response y by defining a ~ U{a,a) and b ^ U{b,b): 



and write the polynomial chaos expansion for the uncertain 
arguments and the response: 

p-i 
a^Y.a,^,{^i) where 5i~W(-l,l) (15) 

1=0 

p-i 
^ = E^J^j(^2) where ^^^^(-l,!) (16) 



3=0 
P-1 

y=^ykipk{$,) where P = 

k=0 



{n+p)\ 
n\p\ 



(17) 



For this paper we are only concerned with the Uniform 
distribution, however due to the specificity of this application 
any other probability distribution with compact support can 
be used (eq. Beta). The sensitivity of the method with respect 
to the shape of the probability density function over the 
compact support remains to be studied. 

Here n is the number of uncertain input variables and 
p the order of the polynomial chaos expansion. The basis 
function ip„ is the n-th degree Legendre polynomial and the 
polynomial coefficients of the input variables are given by: 



flo 



fli 



bo = —z— , oi 



a — a 

2 
b-b 



0-2 



, &2 



flp-l 



Op-1 



(18) 
(19) 



The first six multidimensional Legendre polynomials for 
the bivariate case are given by: 



M^) ^M^i)M^2) =1 

^2(0 =^'2(6)^0(6) =^(3e?-i) 

M^) =M^i)M^2) =66 

1 



(20) 
(21) 

(22) 
(23) 
(24) 



^5(0 = ^'2(6)V^i(6) =2(36-1)6 (25) 

We are interested in finding the polynomial coefficients 
yk which characterize the stochastic behavior of the output 
variable. Using the Galerkin projection and the orthogonality 
property of the polynomials one can isolate the coefficients 
yk as shown in Ref. [5]: 



1 



Vk 



< f, i'k > ^ 



fi^k^{^)d$ 



(26) 



y = fia,b) 



(14) 



where (^(^) = nr=i''«(6) is the joint probability density 
function. The integral can be evaluated using sampling or 
quadrature techniques. 

We show that by bringing the polynomial chaos expansion 
to a Bernstein form using the Garloff's method [8], we can 
efficiently find the minimum and the maximum value of the 
compact support thanks to the properties of the Bernstein 
polynomials: the smallest and the largest coefficient bound 
the output of the function modeled. 

To transform our expansion from Legendre polynomial ba- 
sis to Bernstein polynomial basis we expand our Polynomial 



Chaos expansion, Eq.dTTJ) into a simple power series and 
identify the new coefficients: 



y=Yl ail' 

KN 



(27) 



where the muhi-index I = (ii,...,j„) £ N" and N = 
(ni, . . . , Un) G N" is the multi-index of maximum degrees; 
thus the maximum degree of S^^. is given by rifc. Here, we 
denote ^^ ~ Ci ■ ■ ■ ■ ■ Cn ^^'^ ^he inequality I < J implies 

«1 < jl,---,«7i < in- 

B. Garloff's method to calculate the Bernstein coefficients 

We are interested in the transformation of the power series 
in Eq.dZTli into its Bernstein form: 



y 



E/3iBiN(0 



(28) 



KN 



where Bj^(|) is the Ith Bernstein polynomial of degree N 
on the general box G = [|,|]. In our bi-variate case G = 
[-1,1] X [-1,1] sinceei,6-Z^(-l,l). 



Br(0 = sir(a)--.--si:"(c«) 



(29) 



The univariate Bernstein polynomial i?]? {£,) on the general 
interval [£,,£] is given by: 



BUO 



(e-e)'-(e-e) 



n — k 



The Bernstein coefficients /3i are given by: 

J<I<N \j) 



(30) 



(31) 



where we write (j) = C^) • . . . • (*."). 

The scaled coefficients di are obtained as described in 
Ref.[l] from the ai coefficients in Eq.dZTJi and the box G: 



ai 



ai 



5i(|-0' 



I<J<N 



E i «i 



^J I 



(32) 
(33) 



Example: Consider the following Polynomial Chaos ex- 
pansion given by n = 2,p = 3 and P = 10: 

y^5M0 + M^) + M^) + M^) 

Transforming this expansion into a simple power series 
we obtain the folowing polynomial: 

y = 5 + a + 6 + 66 

where aoo = 5,aio = l,aoi ~ 1. and an ~ 1. Here the 
multi-index of maximum degree is N = (1, 1). 

The following intermediate coefficients are obtain from 
Eq.(l33]l in order to perform the scaling operation: aoo ~ 
4, aoi ~ 0, aio ~ 0, and an = 1. The final set of power- 
coefficients is given by Eq.(l32]i: aoo = 4,aoi = 0, aio = 0, 
and dii =4. 



Finally, the Bernstein coefficients are obtain using Eq. dsTT l: 

/?oo = 4, /3oi = 4, /3io = 4, and (3ii = 8, and the Bernstein 
basis is given by: 

B?,J =i^(l-6)(l-6) Bl\ =1(1-^0(6 + 1) 

^11 =i^(6 + i)(i-6) B}} =l(a + i)(6 + i) 

C. Bounding the range of polynomials 



Given the Bernstein expansion in Eq.(l28]l, the range en- 
closing property [9] gives a bound on the polynomial in terms 
of the Bernstein coefficients: 



min/3i < 2/(1) < max /3j V^eG 



K,l] 



(34) 



Provided that the initial box is small enough, the range 
provided by the Bernstein form is exact. Compared with 
other forms in estimating the range, it is experimentally 
shown in Ref. [25] that the Bernstein form provides the 
smallest average overestimation error in the univariate case. 
For the previous example the range of y is bounded by [4, 8]. 

Tighter bounds can be obtained by subdivision of the 
initial box and choosing the minimum and the maximum 
of all the Bernstein coefficients corresponding to each sub- 
box. An efficient algorithm for range computation that in- 
corporates a number of features such as subdivision, cut-off 
test, simplified vertex test, monotonicity test and others is 
provided in Ref. [20]. 

Therefore, getting back to our problem in mapping DS 
structures on closed intervals, Eq.(fT2]i and Eq.(fT3b, the output 
interval or the focal element Y — [y, y] is given by: 



V = min 01 and 

- kn'^ 



y 



niax/3j 

J<N 



(35) 



This methodology is applied to map all the focal elements 
in the initial body of evidence through the nonlinear function. 
Their corresponding masses are obtained using Eq. (IT2] l. This 
way a body of evidence for the response is constructed. 

IV. NUMERICAL RESULTS 

To prove the concept, we have selected an algebraic 
problem from a set of challenges used to investigate the 
propagation of epistemic uncertainty [17]. The presented 
problem has been investigated previously in the literature 
by Oberkampf and Helton [16]. In the present paper we are 
using the exact parameters for the simulation as in Ref. [16]. 

Consider the following mapping: 



y = f{a,b) = {a + br 



(36) 



where the information concerning a and b is provided by the 
following sources and their corresponding bpa: 



Ai 



A2 



([0.6,0.9], 1.0) 



([0.1, 0.5], 0.2) , ([0.5, 1.0], 0.8) 



B2 



[0.3, 0.5], 0.1) , ([0.6, 0.8], 0.9) 



([0.2, 0.4], 0.1) , ([0.4, 0.6], 0.7) , ([0.6, 1.0], 0.2) 



B3 : I ([0.0, 0.2], i) , ([0.2,0.4],!) , ([0.3, 0.5], i) 

Given the above information, we are looking to bound 
the probabihty of the response in the unsafe region when 
y > 1.7. The function and the desired unsafe region are 
shown in Fig|2] 




In Table U given the ith focal element of A and the j\h 
focal element of B, the Box# is given by 3(j — 1) + i. The 
numbers in bold indicate that a smaller lower bound or a 
larger upper bound has been found by the genetic algorithm. 
In this particular example we have overestimated most of the 
lower bounds and we have provided no underestimation for 
the upper bounds. 

The focal elements from the DS structure in Table U 
are also graphically presented in Fig[3] The wider boxes 
represent the bounds found by the proposed approach while 
the narrow boxes depict the bounds returned by the genetic 
algorithm, and the lines show the range computed using 
interval arithmetics. 
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Induced DS 


STRUCTURE FOR THE RESPONSE y 




Box# 


y. 
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ruf 


Box# 


y 


y 


m.f 


1 


0.687 


0.909 


0.011 


10 


0.890 


1.061 


0.023 


2 


0.721 


1.222 


0.044 


11 


0.967 


1.630 


0.093 


3 


0.741 


1.097 


0.056 


12 


1.007 


1.450 


0.117 


4 


0.804 


0.961 


0.014 


13 


0.953 


1.152 


0.030 


5 


0.853 


1.426 


0.058 


14 


1.069 


1.834 


0.120 


6 


0.880 


1.275 


0.072 


15 


1.123 


1.623 


0.150 


7 


0.850 


1.012 


0.014 


16 


0.952 


1.236 


0.007 


8 


0.912 


1.528 


0.058 


17 


1.068 


2.039 


0.027 


9 


0.945 


1.363 


0.072 


18 


1.122 


1.794 


0.033 



Fig. 2. The function and the unsafe region 

The bpa for a and h is obtained by aggregating the 
information from the first two sources and the last three 
sources respectively using the mixing rule in Eq.® under 
the equal reliability assumption. Thus the new DS structures 
obtained are given by: 



A 



([0.1, 0.5], 0.1) , ([0.5, 1.0], 0.4) , ([0.6, 0.9], 0.5) 



B : <i ([0.0, 0.2], 0.111) , ([0.2, 0.4], 0.144) 
([0.3, 0.5], 0.144) , ([0.4, 0.6], 0.233) 
([0.6, 0.8], 0.3) , ([0.6, 1.0], 0.067) 



The focal elements of the DS structure y are obtained by 
propagating the product space A x B through the nonlinear 
function in Eq.([36]l, and the basic probability assignment is 
obtained using Eq.(fT2b. Thus the induced body of evidence 
for the response y is given in Table U 

The following parameters have been used in obtaining 
the polynomial chaos expansion of the response and the 
afferent bounds: n = 2,p = 5, and P = 19 such that 
the total degree of the polynomial is no greater than 5. 
The integrals in Eq.(l26b have been numerically evaluated 
using Gauss-Legendre quadrature rule with 20 points in 
each direction, and the Bernstein coefficients have been 
obtain using 1 1 subdivisions in each direction. The obtained 
intervals are compared with the intervals given by interval 
analysis (INTLAB [21]), and with the reference intervals 
provided by a genetic algorithm (GA). 




Pl(y>1.7) 



Interval PC-Bernstein 
Reference Interval GA 
Interval I A 



1 2 3 4 5 6 7 8 9 101112131415161718 
Box# 

Fig. 3. Focal elements from the induced DS structure 

To compute the lowest and the highest probability of y > 
1.7, we use the Fig|3]to sum over all the bpa's of the intervals 
that are properly included in the unsafe region for the lower 
bound and for the upper bound, sum over all the bpa's of 
the intervals that intersect the unsafe region. No intervals 
are properly included in the unsafe region, thus the lower 
bound is 0.0. Three intervals are found to intersect the unsafe 
region, boxes: 14, 17 and 18. Summing over their bpa's we 
found the upper bound to be 0.18. All three methods give 
0.0 < Proh{y > 1.7) < 0.18, which is in agreement with 
the result pubUshed by Oberkampf in Ref.[16]. 



Both the CCBF and the CCPF given by Eq.(|9ll-([T0ll are 
plotted in Fig|4]along with the marking for the unsafe region. 
A reason of concern is the overestimation of the lower bound, 
due to the finite polynomial chaos expansion, which in this 
example may provide a larger lower bound for the probability 
of failure. Observe the gross bounds provided by the interval 
arithmetics due to dependency effect. 



Bel : PC-Bernstein 

PI : PC-BernsteIn 

Bel : GA 

PI :GA 

Bel : Interval Analysis 

PI : Interval Analysis 




Fig. 4. CCBF and CCPF from the induced DS structure 



V. CONCLUSIONS 

A new approach in approximating the interval propagation 
for epistemic uncertainty quantification has been presented. 
The input variables are represented as a polynomial expan- 
sion of random variables on compact support, and by apply- 
ing the Galerkin projection in a non-intrusive way we find 
the response of the system also as a polynomial expansion, 
here in the Legendre basis. It exploits the efficient mapping 
of random variables using polynomial chaos expansion from 
which it extracts only the bounds of its compact support. 

We further propose to transform the output polynomial 
chaos expansion from the Legendre basis to the Bernstein 
basis, and use the range enclosure property of Bernstein 
forms to efficiently extract the bounds of the range of 
the output. The method does not suffer of the dependency 
problem as in the interval arithmetics, however much work 
remains to be done in studying the accuracy of interval 
propagation using polynomial chaos. The proposed approach 
is applied on a challenge problem previously investigated by 
Oberkampf et al. to propagate epistemic uncertainty and the 
numerical results obtained provide a basis for optimism. 
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